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Abstract 

Background: Models of biochemical systems are typically complex, which may complicate the discovery of 
cardinal biochemical principles. It is therefore important to single out the parts of a model that are essential for the 
function of the system, so that the remaining non-essential parts can be eliminated. However, each component of 
a mechanistic model has a clear biochemical interpretation, and it is desirable to conserve as much of this 
interpretability as possible in the reduction process. Furthermore, it is of great advantage if we can translate 
predictions from the reduced model to the original model. 

Results: In this paper we present a novel method for model reduction that generates reduced models with a clear 
biochemical interpretation. Unlike conventional methods for model reduction our method enables the mapping of 
predictions by the reduced model to the corresponding detailed predictions by the original model. The method is 
based on proper lumping of state variables interacting on short time scales and on the computation of fraction 
parameters, which serve as the link between the reduced model and the original model. We illustrate the 
advantages of the proposed method by applying it to two biochemical models. The first model is of modest size 
and is commonly occurring as a part of larger models. The second model describes glucose transport across the 
cell membrane in baker's yeast. Both models can be significantly reduced with the proposed method, at the same 
time as the interpretability is conserved. 

Conclusions: We introduce a novel method for reduction of biochemical models that is compatible with the 
concept of zooming. Zooming allows the modeler to work on different levels of model granularity, and enables a 
direct interpretation of how modifications to the model on one level affect the model on other levels in the 
hierarchy. The method extends the applicability of the method that was previously developed for zooming of 
linear biochemical models to nonlinear models. 



Background 

One of the main reasons for the rapid growth of the 
field of systems biology is that it makes extensive use of 
mathematical modeling [1-3]. This allows for a better 
handling of high complexity, which is an inherent prop- 
erty of all living systems. Using modeling, complex 
hypotheses can be formulated and tested in a more sys- 
tematic manner than is possible using only biochemical 
reasoning [4-6]. However, even if one can obtain a 
detailed model of the system with a high predictive 
power, the model in itself does not automatically lead to 
a full understanding of the underlying biochemistry. 
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One should for instance analyze the model to single out 
its essence, i.e., to identify those parts of the model that 
can be eliminated, while still preserving the model's cru- 
cial behavior. This latter task is referred to as model 
reduction, and it is the topic of this paper. There is an 
extensive literature available on the topic of model 
reduction. However, most of these studies have been 
done outside the field of systems biology, and since sys- 
tems biology brings about new types of challenges, 
reduction of biochemical models is still in its early 
stages. Traditional engineering approaches like balanced 
truncation have focused on preserving the input-output 
profile in an optimal manner, both for linear [7-10], and 
for nonlinear [11] systems. However, these methods are 
not suitable for systems biology, because the reduced 
model has no natural interpretation in itself 



O© 201 1 Sunnaker et al; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative 
BiolVlGCl Ccntrsl Commons Attribution License (http://creativecommons.Org/licenses/by/2.0), which permits unrestricted use, distribution, and 
reproduction in any medium, provided the original work is properly cited. 



Sunnaker et al. BMC Systems Biology 201 1, 5:140 
http://www.biomedcentral.eom/1 752-0509/5/1 40 



Page 2 of 21 



(nevertheless, some special cases where this problem can 
be circumvented have been identified [12,13]). This lack 
of interpretation is a problem because systems biology 
models are usually developed to help characterizing the 
dominating parts and structure of the system, and not 
only to obtain a black-box predictor. Methods have 
therefore been developed with traditional chemical 
approaches that are more centered on reducing the 
internal dynamics of the system. These methods are 
typically based on a sensitivity analysis [14-17], on time- 
scale separation [18-21], or on the lumping of state vari- 
ables [22-26] (see [20] for a general review on model 
reduction). The perhaps most widely used method is 
lumping. Two of the main reasons for this are that an 
effective lumping scheme can be identified from basic 
properties of the model (e.g., the stoichiometry), and 
that lumped state variables are formed as easily interpre- 
table pools of state variables in the original model. How- 
ever, lumping does normally not come with the 
possibility of back-translation from the lumped state 
variables in the reduced model to the original state vari- 
ables. In [27] we provided such relations. This means 
that we can take the result from a simulation of a 
reduced model, and without performing a new simula- 
tion, directly compute the corresponding trajectories of 
the desired original state variables. Because of this back- 
translation possibility, we refer to the resulting two 
models as two degrees of zooming of the same model. 
Nevertheless, like in other recent model reduction 
papers in systems biology [28-32], the results in [27] 
were mainly developed with linear systems in mind. Lin- 
ear systems virtually only appear in the cases of mono- 
molecular reaction networks and for models describing 
the probabilistic evolution of a single protein complex 
[27,33]. However, already in [27] we proposed that 
zooming may in principle also be applicable to nonlinear 
models, but we did not derive formulae for back-transla- 
tion. Note that a majority of the currently available sys- 
tems biology models are in fact nonlinear. 

With the method introduced in this paper, we provide 
the extension of the previously proposed method in [27] 
to nonlinear models. We show that new challenges arise 
due to the nonlinearities, but also how these challenges 
can be overcome, for instance with a wise choice of 
state variables in the reduced model. The method is 
demonstrated by application to two closed models of 
metabolic systems. 

Methods 

In this paper we present a more general version of the 
method that was introduced in [27], which is applicable 
to nonlinear models. We start with some basic definitions 
and key observations that are illustrated on a small exam- 
ple model, before we turn to the details of the method. 



Basic Definitions and Assumptions 

The method is developed for models of biochemical 
reaction systems on state space form that are based on 
nonlinear ordinary differential equations (ODEs) 

x = f(x,p,u,t), (1) 

y = h(x,p,u,t), (2) 

where t denotes time; the dot over x in Eq. (1) denotes 
derivative w.r.t. to time; the state vector x e R w ; the 
parameters p e R p ; the inputs u e R m ; the outputs y e 
R l ; and /and h are in general nonlinear functions. The 
state vector, whose individual elements are referred to 
as state variables, typically represents amounts or con- 
centrations of chemical species, and the parameters 
commonly represent kinetic constants, initial conditions, 
or scaling factors. In this paper we are primarily inter- 
ested in a comparison of the state variables between 
models (the original model and the reduced model), 
which means that the form of the nonlinear function in 
Eq. (2) is irrelevant for the application of our method. 
The right-hand side of Eq. (1) can be expressed as the 
stoichiometric matrix S e R nxq times a vector of reac- 
tion rates r = r(x, p, u, t); re R q 

x = Sr (x,p, u, t). 

The existence of separate time-scales are commonly 
utilized for reduction of biochemical models (e.g., by 
reduction of mass action kinetics to Michaelis-Menten 
kinetics). The typical approach is to investigate if subsets 
of the state variables are in steady state or in quasi- 
steady state (QSS). If state variable x t is in steady-state 
for t > 0 it holds by definition that 

Xi(t) = 0, (3) 

which implies that 

^i(0 = ^i(0), 

which efficiently removes the state variable from the 
model, since it can be substituted for constant. If on the 
other hand the state variable x t is in QSS, there are 
terms on the right-hand side of the ODE that are much 
larger than the negligible term on the left-hand side. 
The approximation 

fi(x,p,u,t) « 0. (4) 

is then commonly used to reduce the model. We refer 
to a state as fast in the time interval < T 0 < t < T 1 if Eq. 
(4) is valid in this time interval, and holds for the class 
of all considered inputs to the system. Note that 7\ = °° 
in the case that the systems remains in QSS, which may 
for example not be the case for models with switches 
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(where, e.g., the values of a subset of the state variables 
may change when a certain condition is fulfilled) [34,35]. 

Note that Eq. (3) (steady state) necessitates that Eq. 
(4) (quasi-steady state) is fulfilled, but not vice versa. 
Although QSS implies that (some of) the terms of the 
right-hand side of the ODE are large and leaves the left- 
hand side (derivative term) negligible, the derivative 
term may still be large enough for the state variable in 
QSS to change considerably during the time-span of a 
simulation; the key is that these changes mainly occur 
on a slow manifold. 

Zooming of Linear Models 

The concept of zooming was introduced in [27], and a 
method was presented that is applicable to linear time- 
invariant (LTI) models, which on state space form reads: 



Mo 



x = Ax + Bu 
y = Cx + Du, 



where A e R n x w , B e R KXW ,Ce R l x n , and D e R l 
x m . The method is based on the existence of at least 
one subset of state variables in M 0 for which the inter- 
nal dynamics is very fast with respect to the current 
time-scale of interest. An algorithm for automatic 
reduction of linear models that is based on the detection 
of such subsets, which are referred to as fast clusters, is 
presented in [27]. If the w state variables of a fast cluster 
are replaced by a single state variable x' lL = J27=i x k> we 
obtain a reduced version of the original model 



Mr 



x r = A r x r + B r u 
y = C r x r + Du, 



where x r e R {n ' w + l \ A r e R (n ' w + 1} x {n ' w + l \ B r 

e R (« - w + 1) x m an(J ^ e R / x (* - w + 1) 4 

The fraction parameters, which are typically computed 
from QSS assumptions and mass conservation relations, 
take the form 



(5) 



The fraction parameters are used for back-translation 
of the lumped state variable to the original state vari- 
ables. Note that the fraction parameters are functions of 
the model parameters only, and therefore time-invariant; 
as we will see, these fraction parameter properties do in 
general not hold for nonlinear models. By comparing 
the reactions of the original and reduced models, we see 
that 



(6) 



where k r - L is the rate parameter in the reaction from 
state variable x\ L to state variable Xj in the reduced 
model, and kj t is the rate parameter in the reaction from 
state variable x t to state variable Xj in the original model. 

Finally, note that Eqs. (5) and (6) provide a link 
between Ai 0 an d M r > which constitute two different 
levels of granularity. It is this link between the models 
that make us consider them as two different degrees of 
zooming, and the primary goal of this paper is to estab- 
lish such a link also for nonlinear models. 

Extension to Nonlinear Models - Initial Observations 

We will now present some key observations that are 
used in the derivation of the method for zooming of 
nonlinear models. 

First observe that the mass, which corresponds to a 
weighted (w.r.t. the molecular weight) sum of the state 
variables, of a closed (no exchange of matter with the 
surroundings) nonlinear model is conserved. However, 
the total number of molecules is in general not con- 
served in such a model as it is for linear models. This is 
for example due to the formation and dissociation of 
complexes, which alters the total number of molecules 
in the system. For instance, the binding of A to B 
reduces the number of molecules, as the product AB 
only counts as one molecule; binding reactions cannot 
occur in linear models. A second, and related, observa- 
tion is that another type of conservations appears in 
nonlinear models; conserved moieties. A moiety is a 
specific functional part of a molecule, and the weighted 
sum of the number of molecules that contain this func- 
tional part is constant in a closed system. The presence 
of such a conserved moiety is equivalent to the exis- 
tence of a row vector m e N n for which mS = 0, which 
also implies that 



mx = mSr [x, p) = 0. 



(7) 



If we let the rank of S be denoted by n n the number 
of linearly independent vectors for which Eq. (7) holds 
is equal to n - n r , which implies the existence of a 
matrix M 



M S = 0, 



(8) 



i=i 



where M e [N n-nrXn - 

Let us now make some remarks regarding fast state 
variables in a nonlinear model. Let Xf e R n f be the vec- 
tor of all fast state variables in T 0 < t <7\. For simplifi- 
cation we will assume that there are no inputs to the 
system, although it would in principle be possible to 
incorporate inputs in the following discussion. The 
right-hand side of the ODEs for these fast state vari- 
ables, if there are no inputs, can be separated into two 
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parts. The first part contains reactions between fast state 
variables that are significant for the fast dynamics; Vf (xp 
p), and the second part contains all other reactions, r s (x, 
p), i.e., 



kf = S s r s (x,p) + SfTf[xf,p), 



(9) 



where Sf and S s are the corresponding stoichiometric 
matrices. Let us now consider the fast stoichiometric 
matrix, Sp and especially the conserved moieties that are 
implied by Sf, Since these moieties are only (approxi- 
mately) conserved on a fast enough time-scale, we refer 
to such moiety conservations as apparent conservations. 
Let Mf be a matrix with a linearly independent rows 
such that 



MfSf = 0, 



(10) 



where Mf e h\ axn f. Each row of this matrix thus 
implies an apparent conserved moiety in the system. Let 
the sums of state variables that correspond to apparent 
moiety conservations (i.e., lumps of state variables) be 
denoted by /, so that 

l = M f Xf. (11) 

If we differentiate / with respect to time, we get 



/ = Mfkf = Mf(S s r s (x,p) + SfTf(xf,p)) 
^M f S s r s (x, p), 



(12) 



where Eqs. (9) and (10) were used. 

It is interesting to note that the matrix Mf is not unique, 
but that in fact any matrix M = NMf can be used for 
lumping, where N e R axa is non-singular. This observa- 
tion allows us to choose a matrix £4 for which a maximal 
number of rows in M S s vanish, which results in the great- 
est possible reduction in the number of state variables. 
Finally note that Eq. (4), in the absence of inputs, gives 



ff(x, p, u t , t) = ff(x, p, 0, t) = 

= S s r s (x,p) +SfTf(xf,p) 
^S f r f {x f ,p)^0 



(13) 



since the term S/y {xp p) dominates the term S s r s (x, p). 
Note that Eq. (4) and consequently Eq. (13) only hold in 
T 0 < t <T lf since the state is only known to be fast in 
this time span. 

Eq. (12) defines the ODEs of the reduced model (the 
lumped state variables), and Eqs. (11) and (13) can in 
principle be used to calculate back-translation formulae, 
as is demonstrated with the small example model in the 
next section. However, as we shall see, this approach 
requires the explicit algebraic solution to a system of 



nonlinear equations, which is typically an infeasible task. 
Furthermore, there is not a clear one-to-one mapping 
between the state variables of the original and reduced 
models as in the case of proper lumping [27]. 

A Small Example Model 

We will now present a small example model, with three 
fast state variables, which is reduced with the approach 
discussed above. An alternative approach is then 
demonstrated with the advantage that it scales better to 
larger models. 

Consider the reversible formation of a complex C 
from a substrate A and an enzyme B 



A-\- B 



consisting of the fast state variables Xf = {A B C) r , 
where the bullets (♦) represent the slow state variables 
surrounding the three fast state variables in the model. 
The ODEs for the fast state variables take the form 



Xf ■ 




(feiAB - fe-iC) + S s r s {x,p). 



(14) 



The three state variables in the model constitute a fast 
cluster with two apparent conserved moieties, which 
may be represented by the following relations 




: MfXf, 



(15) 



where the lumped state variables L x and L 2 are intro- 
duced. Note that Eq. (12) defines the dynamics of the 
lumped state variables. 

The distribution of mass among the fast state variables 
is given by Eq. (15) and by applying Eq. (13) to (14), 
which results in an equation system with the three fast 
state variables as unknowns 



C = Li 



B + C = L 2 . 



(16) 



(17) 



(18) 



Analytic expressions for the fast state variables A, B, 
and C are given by the non-negative solution to Eqs. 

(lexis) 

A « h U _ L2 _ Kl + y (Ll +L 2 +K 1 ) 2 -4L 1 L 2 )i;i9) 
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I(_ Ll + L2 _ Kl + y /(L 1 +L 2 + K 1 ) 2 -4L 1 L 2 ), (20) 



C « — (Li + L 2 + + y (Li + L 2 + Ki) 2 - 4LiL 2 ), (21) 



where K"i 



fel ' 



We can now employ Eq. (12) to solve the ODEs for 
the lumped state variables L x and L 2 , and use Eqs. (19)- 
(21) as back-translation formulae to compute the trajec- 
tories of the original state variables A, B, and C. How- 
ever, note that for even slightly larger clusters of fast 
species than the one discussed here it would not be pos- 
sible to calculate algebraic expressions of the original 
state variables with this approach, since it builds on the 
explicit solution of a system of nonlinear equations, 
which quickly becomes infeasible with growing problem 
size. 

Alternatively, we can take an approach to the problem 
that is inspired by the method for linear systems in [27]. 
The first step is to express Eqs. (16) and (17) as a linear 
system w.r.t. the state variables A and C 



kiB -fe_i 
1 1 



The solution to Eq. (22) w.r.t. A and C is 



A= y]a{B,P)U 



C = ri c {B,p)Li 



fe_ 



Ki 



B + K 1 



B 



B + K\ 



(22) 



(23) 



(24) 



where K\ = — -, and the fraction parameters T] A (B t p) 

and r\c (B, P) are defined in Eqs. (23) and (24), respec- 
tively. The ODE for L x is defined by Eq. (12), and the 
ODE for B can be derived by differentiation of L 2 in Eq. 
(18), which gives that 



dB 
dt 



(1 + 



r)" 1 



B 



(B + Ki) 2J V B + K 
KiLi / B \ 

MfS s r s (x,k), 



t 
\dtJ 



(25) 



The reduced model consists of the two state variables 
Li and B (note that L 2 does not appear in the reduced 
model), and the dynamics is described by Eqs. (12) and 
(25), respectively. Note that the state variables A and C 
can be back-translated from the reduced model with 



Eqs. (23) and (24). This approach is a bit more intricate 
than the first, but comes with the advantage that we do 
not need to solve a system of nonlinear equations. 

A Method for Zooming of Nonlinear Models 

We will now step-by-step present a method that can 
be used to construct zoomable nonlinear biochemical 
models. This involves two sub-goals: i) to identify a 
reduced model that shares important characteristics 
with the original model, ii) to derive back-translation 
formulae that can be used to compute the original 
state variables and parameters from the reduced 
model. 

In an initialization step of the method for a model 
M 0 we first formulate mathematical equations for all 
conservation relations Eq. (8), state variables in steady- 
state Eq. (3), and quasi-steady state assumptions Eq. (4). 
If additional properties of the system are known, we 
also formulate the corresponding equations. 

Step 1 

The first step of the method is to identify the apparent 
conservation relations in the model. 

Definition 1: Let Sf be the stoichiometric matrix for 
the reactions jy (Xf p) as defined in Eq. (9). Each subset 
of state variables for which the corresponding rows of Sf 
are linearly dependent constitutes an apparent conserva- 
tion relation. Hence the apparent conservation relations 
lie in the left null space of Sf and the dimension of this 
space is n - rank(^). 

Note that the apparent conservation relations are 
defined in Eq. (11). It is trivial to identify the set of all 
linearly dependent rows of Sf with a mathematical com- 
puting software (e.g., SBtoolbox for Matlab [36]). 

Step 2 

The second step of the method is to define the state 
variables of the reduced model, which we refer to as 
modified lumped state variables. 

Definition 2: Let x be a lumped state variable corre- 
sponding to a subset of the state variables in an appar- 
ent conservation relation. Then x is a modified lumped 
state variable if the lumping scheme with respect to the 
state variables of the original model is proper. 

Note that the original state variables have a clear 
interpretation in the reduced model (i.e., that the 
lumped variables form disjoint sets) if the lumping 
scheme is proper, i.e., 



l m — bA m Xj 



(26) 



where M m is a x rif matrix with elements equal to 0 or 
1 and column sums equal to 1, and l m denotes the mod- 
ified lumped state variables. We typically have a large 
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freedom in the choice of M m . The number of state vari- 
ables is maximally reduced if all exact conservation rela- 
tions in the model are retained as modified lumped 
state variables (and replaced by constants). 

Step 3 

The third step of the method is to derive fraction para- 
meters, which constitute the link between the reduced 
model and the original model Let the original state vari- 
ables that constitute the /c:th modified lumped state vari- 
able l m be denoted by x mk , so that 



i=l 



mki' 



(27) 



A number of n m equations that are linear w.r.t. x mk , 
and linearly independent, are required to calculate frac- 
tion parameters. The existance of n m such equations 
results in an equation system 



h{lm,p) = A(l m/ p)x n 



(28) 



where both A{l m ,p) e R n ™ XUm and bk{l m ,p) e R Hm are 
known, although some of the equations may in general 
be approximate (e.g., QSS). The matrix A(l m , p) is inver- 
tible since the equations are linearly independent, and 
we have that 



Xm k = A 1 (lm,p)h(!m,p)' 

The fraction parameters can then be calculated 



r)m hi {lm,p) 



i=l x m u 



(29) 



(30) 



where we used Eq. (27) in the last step. 

A modified lumped state variable for which an insuf- 
ficient number of linear and linearly independent 
equations are available may still be used in the 
reduced model. However, the back-translation of the 
modified lumped state variable to the original state 
variables is then not possible, and step 3 of the 
method is ignored. 

Step 4 

The fourth step of the method is to derive the rate of 
change of the modified lumped state variables. Theo- 
rem: The dynamics of the modified lumped state vari- 
ables is given by 



\ m = (I + J(l mf p))- 1 'l = 

= {I + ]{l mi p)T l M f S s r s {x l p) l 



(31) 



where S s and r s (x, p) were defined in Eq. (9), the 
matrix Mf is defined in Eq. (11), and 



hj{lmrp) = Y,^ M h- M rn lk ) 



dgWm,p) 

dim. 



(32) 



where the matrix M m is defined in Eq. (26), and 
gi{lmrp) = Xf. is introduced to simplify the notation. 
Proof: 

First subtract Eq. (26) from Eq. (11) 

/ =l m + {Mf - M m )x f = l m + {Mf - M m )g(l m , p), 
and differentiate / with respect to time, which gives 

i='lm+ J{lm,P)lm = U + J(lm,p))lm, ( 3 3) 

where / is the identity matrix and /(/ m , p) is the Jaco- 
bian of {Mf - M m )g(l m , p) with respect to l m . The ele- 
ment Jij of /(/ m , p) is given by 



Jit l m>P) = J2( M f*- Mm ^ 



dgk{lm,p) 

dim. 



From Eq. (33) it is straight-forward to derive / m , 
which takes the form 

/m = (/ + /(/m,p))" 1 / = 

= {I + Klm,p)r 1 M f S s r s [x,p), 

where Eq. (12) was used in the last step. D 
The matrix / + J(l m ; p) is symbolically invertible, but 
may in general contain singularities for particular com- 
binations of parameters values and state variable values. 
However, the matrix is always invertible for the models 
discussed in this paper, since the corresponding deter- 
minants are strictly positive. 

Step 5 

The final step of the method is to back-translate the 
modified lumped state variables to the original state 
variables with the fraction parameters derived in step 3. 
This allows a comparison between the predictions by 
the reduced model to those of the original model. 

The implementation of the method is straight-forward, 
and we have used Matlab (R2008b) together with the 
SBtoolbox [36] as computing software for the models in 
this paper. 

Results 

We will now demonstrate the method through applica- 
tion to two example models. 

Enzyme Kinetics Model 

The model below describes the process of conversion of 
a substrate, S, into a product, P, which is catalyzed by 
an enzyme, E. 
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s ... K ' 



P + E 



Note that the complexes C s and C p are formed by S 
bound to E, and P bound to £, respectively. This model 
is frequently occurring as part of larger models of biolo- 
gical systems, although the reaction from C s to C P is 
sometimes neglected, or reversible. The ODEs for the 
model are listed in Appendix A.l, where the three reac- 
tions are defined as: r x - k\SE - k^C S) r 2 - k 2 C s , and r 3 
= k 3 C P - Ic 3 PE. 

The reaction terms ri(x, p) and r 3 (x, p) are assumed to 
be dominating, and the reaction term r 2 (x, p) to be 
insignificant in the ODEs. This results in that all state 
variables are in QSS, which gives 



hSE-k-^s - 0, 



k 3 PE-k- 3 C P « 0. 



(34) 



(35) 



We denote the sum of the state variables containing 
the enzyme by 



L E = E + Cs + Cp, 



(36) 



which is constant since the total amount of the 
enzyme E is conserved in the system. 
Reduction of the Enzyme Kinetics Model 

The first step of the method is to identify the apparent 
conservation relations from the matrix Sf, Since r 2 (x, p) 
is dominated by r x [x, p) and r 3 (x, p) the model ODEs 
can be written on the form of Eq. (9) 



E 
P 
Cs 
\CpJ 



■- S s r s {x,p) +S f r f {xf r p) ■■ 



( 0 ^ 




f-1 0 \ 


0 




-1 1 


0 


r 2 + 
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A basis of the left null space of Sf is given by the row 
vectors of Mp which is defined by 
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(37) 



where L s and L P are apparent conservation relations 
and LE is an exact conservation relation. 

The second step is to define the modified lumped 
state variables on the form of Eq. (26). The number of 
state variables is maximally reduced if L E is retained as a 
state variable in the reduced model (i.e., since L E , unlike 
L s and L P) can be replaced by a constant). The vector of 
modified lumped state variables then is defined 







( s \ 




/l 0 0 0 0\ 


E 




=00100 
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L E J 


\0 1 0 1 I J 


C s 









= MfXf, 

In the third step of the method we calculate fraction 
parameters for the modified lumped state variable L E . 
There are five equations (Eqs. (34)-(35) and (37)) that 
are linear w.r.t. the state variables £, C s , and C P) which 
are lumped into the state variable L E . Note that only n m 
= 3 equations are required to derive fraction parameters, 
and we use Eqs. (34) -(36) to formulate an equation sys- 
tem as in Eq. (28) with the solution 



E = r] E L E 



1 



1 +M 1 S + M 3 P 



Cs = rjc s L E * 



Cp = r] Cp L E 



1 +MiS + M 3 P 



M 3 P 
1 +A4iS + A4 3 P 



Le =g4{lmrp)r 



L E =g5{lm f p) f 



(38) 



(39) 



(40) 



where M\ = ^ and M 3 



The two remaining 

modified lumped state variables correspond to S and P 
in the original model, so we define that 
Um = S = gi{l m ,p) and \ mi = P = g 3 {l mi p) . 

In the fourth step we derive the rate of change of the 
modified lumped state variables. Eq. (32) gives that 
L E = 0> which is replaced by a constant, and 



-k 2 MiSL E {{\ + MiS + M 3 Pf + M 3 L E ) 
0(S, P, L E ,Mi,Ms) ' 



fe 2 A4!SL £ ((l + MiS + M 3 Pf + MiL E ) 



(t>{S,P,L E ,Mi,M 3 ) 



(41) 



(42) 



where 



: MfXf, 



0(S,P,L £ ,A4i,A4 3 ) = ((1 +MiS + M 3 P) 3 +(Mi+M3+MiM 3 (P+S))(l+MiS+M3P)L£+MiM 3 L^). 
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The two ODEs in Eqs (41) -(42) define the dynamics of 
the state variables in the reduced model. We finally note 
that the exact conservation relation for the substrate, L T 
= S + P + C s + Cp, together with Eqs. (39)-(40) can be 
used to reduce the model further to a single state. 

In the fifth step of the method we use the fraction 
parameters, defined in Eqs. (38)-(40), to back-translate 
the modified lumped state variables to the state variables 
of the original model. A comparison between predic- 
tions of the original state variables, from simulations of 
the original model and the reduced model, is presented 
in Figure 1. Implementations of the original model 
(Additional file 1), the reduced model (Additional file 2), 
and a script for simulation with SBtoolbox2 for 
MATLAB [36] (Additional file 3), are available in Addi- 
tional files. 

The only assumption that was used in the derivation 
of the reduced model is that the reaction terms r^x, p) 
and r 3 (x, p) dominate the reaction term r 2 (x, p), which 



results in that all state variables are in QSS. To assess 
the impact of these assumptions on the reduced model 
we compute the relative difference between the state 
variables in the original and in the reduced model 



Si(t) = 



4(0 ' 



i = 1, 



where x- is state variable i in the original model, x\ is 
the corresponding back-translated state variable in the 
reduced model, and \x\ denotes the elementwise abso- 
lute values of x. The maximal mean and infinity norm 
of 8i(t) in Eq. (43) over time is presented in Table 1 for 
parameter values over five orders of magnitude. In gen- 
eral, the reduced model appears to be robust to changes 
in the parameter values, although slightly more sensitive 
to some parameters (e.g., small values of large values 
of ki, or large values of k 2 , which violate the assump- 
tions used in the reduction). However, note that the 
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Figure 1 Small example model. A comparison between the state variables of the original enzyme kinetics model and the backtranslated state 
variables of the reduced version of the same model. 
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Table 1 Robustness of the reduced model for large deviations from the nominal parameter point are presented for 
the enzyme kinetics model, with a sampling frequency of 0.1 (starting from 0.1) time units. 



Param./Factor 10" 2 10" 1 10° 10 1 10 2 



*1 


0.00053/0.0082 


0.0018/0.0071 


0.0010/0.0059 


0.019/0.18 


0.19/1.8 




0.19/1.9 


0.0061/0.079 


0.0010/0.0059 


0.0017/0.0025 


0.00024/0.0028 


K 2 


0.00019/0.0060 


0.00082/0.0055 


0.0010/0.0059 


0.035/0.23 


0.21/0.39 


k 3 


0.0068/0.014 


0.0053/0.010 


0.0010/0.0059 


0.0067/0.044 


0.035/0.29 


K-3 


0.020/0.30 


0.0048/0.032 


0.0010/0.0059 


0.0051/0.0093 


0.0069/0.015 


All 


0.012/1.2 


0.0014/0.062 


0.0010/0.0059 


0.025/0.41 


0.010/0.045 



The nominal parameter values {ki - 1000, /c. q = 2000, k 2 = 1, k 3 = 1000, and /c_ 3 = 3000) are modified by a multiplicative factor, and the maximal (for any state 
variable) time average/infinity norm of the relative difference between the original and the reduced model is presented above. Note that only concentrations 
larger than 10~ 6 are considered in the analysis above, due to potential numerical inaccuracies. 



validity of the QSS assumption may also depend on the 
state variables, for example the total concentration of 
the enzyme. Interestingly, we observed that the reduced 
model can well approximate the original model over 
several order of magnitudes around the nominal enzyme 
concentration (L E = 1). It is well-known that the QSS 
approximation is only valid for sufficiently small enzyme 
concentrations, and as expected the performance of the 
reduced model starts to decrease for immense enzyme 
concentrations. 

Note that all the state variables of the original model 
have a direct biological interpretation also in the 
reduced model, and that Eqs. (38)-(40) can be used to 
back-translate the state variables. The reduced model 
may be depicted 

fc-l fc-3 

■ s • m Li: ^ilc s Li: — — r/, •, /./ _ , /' • /// /./• 

where the fraction parameters specify the distribution 
of the enzyme among the corresponding original state 
variables. 

Glucose Transport in Budding Yeast 

A model for the transport of glucose into a cell of 
baker's yeast (S. cerevisiae), which constitutes the first 
step of glycolysis, is presented in [37]. The inflow of glu- 
cose is modeled as a facilitated diffusion process, in 
which a carrier enzyme is responsible for the transport 
between the inner and outer regions of the cellular 
membrane. It is assumed that glucose 6-phosphate 
(G6P) has an inhibitory role in the glucose transport 
process by binding to the transporter. A graphical repre- 
sentation of the model is shown in Figure 2, and the 
ODEs for the state variables are listed in Appendix A.2. 

In [27] we described how the calculation of fraction 
parameters, based on a set of assumptions, leads to the 
same reaction rates in the reduced model as were 
reported in [37]. The assumptions are that state vari- 
ables participating in reactions for uptake and release of 
glucose and G6P across the cell membrane are in QSS, 
that the transporter is conserved, and that the 



concentrations of the transporter in the inner and outer 
regions of the cellular membrane are constant. 
The assumption that the state variables 

x ac> x Glc> 4-G6P> and 4-G/C-G6P* which participate in 
the uptake and release of G6P and glucose across the 
cell membrane, are in QSS gives that 




Figure 2 Glucose transport model. The original model for glucose 
transport in baker's yeast (S. cerevisiae). This figure was originally 
presented in [27]. 

k ) 
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k 3 X l E _ 



Glc^G6P 



k- 3 x l E _ 



Glc—G6P 



0. 



(46) 



We have the following exact conservation relations in 
the model 



Le ~ x E-G6P + X E-Glc-G6P + X E-Glc + ' ' ' 
X E-Glc + X E + X E' 



Ldc - x c1r + x cl + X ] 



Glc ^ ^Glc ^ ^E-Glc-G6P + X E-Glc + X E-Glc 



LG6P = XG6P + X \ 



X E-G6P 



+ X 



E—Glc—G6P' 



(47) 



(48) 



(49) 



where L E) L Gtc , and L G6P are constant over time. The 
assumption in [37] that the concentrations of the trans- 
porter in the inner and outer regions of the cell mem- 
brane are constant is formulated 



a(x\ 



E-Glc 



■X 



■E-Glc 



) + /H4-4) = o. 



(50) 



It is not clear how equations for back-translation of 
the state variables in the reduced model in [27,37] can 
be derived. The reduced model has three state variables; 
external- and internal glucose and G6P, but only two 
differential equations for the in- and outflow of glucose, 
since the ODE for G6P is replaced by a representative 
function that is inferred from the G6P data. Our 
method does not rely on that such information is avail- 
able, although it would in principle be possible to utilize 
data fitted functions for the state variables. Also note 
that the equations in the reduced model describe the 
total influx and efflux of glucose across the membrane 
[27], which cannot be interpreted w.r.t. the state vari- 
ables of the original model. Other assumptions in [37] 
that complicates a comparison with our method are that 
the efflux of glucose is negligible, that the concentration 
of glucose in the cytosol is negligible, and that the con- 
centrations of the transporter are constant in the inner 
and outer regions of the cell membrane (Eq. (50)). 

It is not possible to generate a reduced model by 
direct substitution of the fraction parameters that were 
derived in [27] into the ODEs of the original model, 
since this would lead to the prediction that the state 
variables are constant (as discussed in [27]). We will 
now instead illustrate how our method can be used to 
derive a reduced and zoomable version of the glucose 
transport model. 

Reduction of the Glucose Transport Model 

Before applying our method to the glucose transport 
model we tried an alternative approach. Eqs. (43)-(46) 
were solved w.r.t. the state variables in QSS, and the 
resulting expressions were then substituted into the 
remaining ODEs. The details of the derivation of the 



reduced model are presented in Appendix A. 3. The 
reduced model does not produce satisfactory predictions 
for any other state variable than x l G6P > which remains 
approximately constant during the simulation. Imple- 
mentations of the original model (Additional file 4), the 
reduced model (Additional file 5), and a script for simu- 
lation with SBtoolbox2 for MATLAB [36] (Additional 
file 6), are available in Additional files. 

Since the first approach turned out to be insufficient 
for reduction of the glucose transport model we 
applied our method to the same model. Following [37], 
we initially assumed constant transporter concentra- 
tions in the inner and outer regions of the cellular 
membrane, as defined by Eq. (50). For the details on 
the derivation of the reduced model we refer to 
Appendix A. 4. The reduced model clearly performs 
better than the model resulting from the first 
approach, but it is still not satisfactory. However, the 
assumption of constant regional concentrations of the 
transporter may not be valid since the transport of glu- 
cose across the cell membrane is a rate limiting step in 
the model, and appears to be important for the state 
variable dynamics. We therefore decided to neglect Eq. 
(50) in the reduction process. 

Implementations of the original model (Additional file 
4), the reduced model (Additional file 7), and a script 
for simulation with SBtoolbox2 for MATLAB [36] 
(Additional file 8), are available in Additional files. In 
the first step of the method we identify the following 
apparent conservation relations 

L G ic 2 =Mx = 

LG6P 

\Le 2 J 

/10000100 0\ 
00000 10 10 
0 10 10 0 10 0 
0 0 1110 0 0 0 
\00 1 1 00 1 0 1/ 

We note that there are two disjoint clusters of fast 
reactions in the model, corresponding to the outer- and 
inner parts of the cell membrane. 

In the second step we define the modified lumped 
state variables. We decide to keep the lumped state 
variables L El and L G6P as modified lumped state vari- 
ables. The choice to keep L G6P leads to the largest pos- 
sible reduction in the number of state variables, since 
the conservation of L G6P is exact (which is not true for 
any other state variable in /). The modified lumped 
state variables are defined true for any other state 



(51) 
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variable in /). The modified lumped state variables are 
defined 



Glc 
i 

Glc 
LG6P 



X 



■ M m x ■■ 



\Le 3 ) 

(I 0000000 0\ 
0 0 0 0 0 10 10 
010000000 
0 0 1110 0 0 0 

\0 00000101/ 



(52) 



where 



X Glc 

X Glc 
X E-G6P 
—GJ.C—G6P 



x E 



G6P 
Glc 



X E-Glc 
X E 



\ 



I 



Note that two of the state variables in the original 
model, x e Gk and x l Glc , are also modified lumped state 
variables. 

In the third step of the method we calculate fraction 
parameters for the modified lumped state variables that 
correspond to more than one of the original state vari- 
ables (i.e., L £l , L G6P , and Le 3 ). All of the modified 
lumped state variables satisfy the requirement that at 
least n m of Eqs. (43)-(49) and Eq. (52) are linear, and 
linearly independent, with respect to the corresponding 
original state variables. Eqs. (43) and (52) form a non- 
linear equation system with the solution 



E-Glc 



^E-Glc 

\K l+ x Gk ) \x Gk 



(53) 



Let us define that #6 = x e-gIc an d #8 = x e- Similarly, 
the fractions of the two carrier state variables in the 
inner regions of the cell to the lumped state variable 
L Es can be computed from Eqs. (44) and (52) (Eq. (28)) 



V E-Glc 



Ve-GIc 



Le 3 



1 



(x ! Gfc + K 2 ) 



K 2 



(54) 



Le 3 



We define that #7 = x l E _ Gk and g 9 = x l E . The fraction 

parameters for the G6P-state variables can be computed 
from Eqs. (45)-(46) and Eq. (28) with the solution 



(55) 



/ X E-G6P ^ 




' ^E-G6P ^ 




1 X E-Glc-G6P 


1 = 1 


Ve-GIc-GGP 


LG6P 


\ X G6P / 




V ^L 6 p / 





K^x E 
K4X l E _ Gk 
K 3 K 4 



Lg 



6P- 



where £ = JC 3 JC 4 + K4 x E _ Gk + Ki x E ♦ We define that 

g 4 = x E -Glc-G6P> # 4 = X E-Glc-G6P> anC * ^ ~ X G6P' F ° r 

the two original state variables that are kept as modified 
lumped state variables in the reduced model we define 
that gi =x e Gk and g 2 =x i Gk . 

The fourth step of the method is to derive rate equa- 
tions for the modified lumped state variables. Since the 
apparent conservations are separated into two disjoint 
clusters of fast state variables, we can treat the model 
for the inner and outer regions of the membrane sepa- 
rately. Let the modified lumped state variables corre- 
sponding to the outer region of the cell membrane be 
denoted by l mi = (x e Gk L £l ) T , and the variables in the 
inner region of the cell membrane by l mi = [x l Gk Le 3 ) t . 
Note that the state variable L G6P can be replaced by a 
constant in the model, since l G6P = 0. The ODEs of the 
modified lumped state variables are derived with Eq. 
(31). In the inner region the ODEs are 



In 



where 



1 +/n Jn 
0 1 



(56) 



LGId 



-a(x\ 



EGlc 



EGlc 
X EGlc 



x EGlc 



)-/H4-4)7' 



(57) 



In the larger outer region of the cell membrane the 
ODEs take the form 



^3 / 



1 + Jn J12 
J 21 1 + J22 



(58) 



h, 
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where 



( Lck 2 \ 
\Le 2 J 



-a{x\ 



EGlc 



EGlc 
' X EGlJ 



X EGlJ 



■ x e E ) 



(59) 



In the fifth step of the method the reduced model, 
which is defined by the four ODEs in Eqs. (56)-(59), is 
simulated. The trajectories of the state variables of the 
reduced model can be back-translated to the original 
state variables with the fraction parameters defined in 
Eqs. (53), (54), and (55). The simulation results are 
shown in Figure 3 and Figure 4. All the state variables 
can be back-translated properly, which shows that the 
model properties that are important for recovery of the 
state variables are retained in the reduction. Implemen- 
tations of the original model (Additional file 4), the 
reduced model (Additional file 9), and a script for simu- 
lation with SBtoolbox2 for MATLAB [36] (Additional 
file 10), are available in Additional files. 

If we use L El instead of L G6P as a modified lumped 
state variable, the reduced model will have the same 



state variables as the reduced model in [27,37] (i.e., 
x Gk' x Gk' anc ^ x G6p) two additional state variables for 
the transporter. This gives a reduced model with five 
state variables, but equally many parameters as in the 
previous case. A comparison between the original model 
and the reduced model, w.r.t. the original state variables, 
is shown in Figure 5 and Figure 6. As can be seen the 
comparison is very good, in fact it is even slightly better 
than for the reduced model with four state variables. 
The details of the derivation of the reduced model are 
presented in Appendix A.5. Implementations of the ori- 
ginal model (Additional file 4), the reduced model 
(Additional file 11), and a script for simulation with 
SBtoolbox2 for MATLAB [36] (Additional file 12), are 
available in Additional files. 

Note that the only assumption used to derive the 
reduced model is that states that are involved in reac- 
tions at the membrane are in QSS. To investigate the 
parameter space region in which the QSS assumptions 
are valid we use the measure defined in Eq. (43). The 
maximal mean and infinity norm of the relative differ- 
ence between the original and the reduced model in Eq. 
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Figure 3 Reduction with our method to four state variables. A comparison between the original glucose transport model and the model 
reduced to four state variables with our method, w.r.t. the state variables of the original model. 
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Figure 4 Reduction with our method to four state variables. A comparison between the original glucose transport model and the model 
reduced to four state variables with our method, w.r.t. the state variables of the original model. 



(43) over time is presented in Table 2. The reduced 
model appears to be relatively robust to changes in the 
parameters, although sensitive to small values of /c 4 and 
to large values of /c 4 . This is mainly due to that a large 
proportion of the transporter E is absorbed in x l E _ G6P , 
which leads to that some of the QSS assumptions are 
invalid. We also observed that relative difference 
between the models is insensitive to the total concentra- 
tion of the transporter for several orders of magnitude 
around the nominal value (L E = 0.01). However, note 
that this observation is specific to the studied model 
and may not be generalizable to other similar biochem- 
ical models. 

Discussion 

In this paper we have presented a novel method for 
reduction of biochemical models that is compatible with 
the concept of zooming. Several methods for reduction 
of biochemical models already exist in the literature. 
However, few of these methods result in biochemically 
interpretable models, and to our knowledge there are no 
nonlinear lumping methods for which the state variables 



and parameters of the reduced model can be back-trans- 
lated (mapped) to the original model. 

The application of the QSS assumption has been a 
commonly used tool in the modeling of biochemical 
networks since the late 1960s, and in chemical kinetics 
for more than 80 years [38]. The validity of the QSS 
approximation is well studied both for specific biochem- 
ical mechanisms [38,39] and for more complex models 
[40,41]. The resulting equations, together with conserva- 
tion relations, are typically used to eliminate some of 
the state variables in the model (e.g., see [40]). However, 
with the examples in this paper we have showed that 
such an approach is not always sufficient, and we pro- 
pose to use proper lumping of state variables in combi- 
nation with back-translation. 

Our method has several important advantages when 
applied to biochemical models. The most important 
advantage is that we end up with reduced models with a 
clear biological interpretation, meaning that each state 
variable of the original model corresponds to a fraction 
of exactly one of the state variables in the reduced 
model. A consequence is that neighboring species in the 
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Figure 5 Reduction with our method to five state variables. A comparison between the original glucose transport model and the model 
reduced to five state variables with our method, w.r.t. the state variables of the original model. 

V J 



original model remain neighbors in the reduced model 
Hence we can consider the original and reduced models 
as two different degrees of zooming; a concept that we 
discussed in some detail in [27] for linear models. 

The work in this paper can be seen as an extension of 
the theory introduced for linear models in [27] to non- 
linear models. The method is based on assumptions 
regarding the dynamics that result in a sufficient number 
of equations that are linear w.r.t. the state variables to be 
back-translated. Such equations are typically a natural 
result of QSS assumptions and conservations relations in 
models based on mass action kinetics [42], and in particu- 
lar in models that involve transporters and enzymes (e.g., 
the models in this paper). However, note that our method 
may also be applicable to models with other types of reac- 
tion kinetics. We also note that if too few linear relations 
are available for calculation of fraction parameters for a 
part of a model, this part can still be reduced and the 
reduced model can be simulated, although we cannot 
back-translate the corresponding modified lumped state 
variables since no fraction parameters are available. How- 
ever, depending on the purpose of the model it may be 
enough to calculate fraction parameters for a subset of the 



state variables in the reduced model. Linearization of the 
model around a steady state operating point may also be a 
feasible approach to calculate fraction parameters with the 
method in [27]. 

The proposed method enables mapping of the state 
variables and parameters of the reduced model to those 
of the original model. In [27] we referred to this mapping 
as back-translation. Back-translation is of great impor- 
tance, since we can directly observe how modifications to 
the reduced model impact the original model. It also 
gives the modeler an opportunity to check whether the 
assumptions underlying the reduction are acceptable. To 
illustrate the power of back-translation we provide plots 
for comparison of simulations of the original and reduced 
models, w.r.t. the original state variables, for the models 
to which the method is applied in this paper. 

Back-translation of state variables typically requires 
the solution of a system of nonlinear equations, which 
often results from the assumption of state variables in 
QSS and conservation relations. Unfortunately, analytic 
solutions to systems of nonlinear equations do in gen- 
eral not exist. An advantage with the proposed method 
is that such solutions are not required, since they are 
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Figure 6 Reduction with our method to five state variables. A comparison between the original glucose transport model and the model 
reduced to five state variables with our method, w.r.t. the state variables of the original model. 

V J 



replaced by computation of the inverse of a matrix for 
each cluster of fast state variables, which is in general a 
more feasible task. 

Our method was applied to a small model with five state 
variables that commonly appears as part of larger 



biochemical models, and to a previously published model 
for the transport of glucose in baker's yeast (S. cerevisiae) 
[37]. The first model was reduced from five to one state 
variable, and from five to three parameters. However, note 
that our focus has been on the reduction of the number of 



Table 2 Robustness of the reduced model for large deviations from the nominal parameter point are presented for 
the glucose transport model, with a sampling frequency of 1 (starting from 1) time units. 



Param./Factor 10' 2 10' 1 10° 10 1 10 2 



/C1 


0.079/0.091 


0.079/0.10 


0.081/0.091 


0.083/0.10 


0.084/0.1 1 




0.085/0.23 


0.084/0.093 


0.081/0.091 


0.078/0.10 


0.078/0.089 




0.16/0.93 


0.080/0.55 


0.081/0.091 


0.063/0.083 


0.054/0.062 




0.054/0.10 


0.062/0.074 


0.081/0.091 


0.097/0.10 


0.10/0.10 




0.081/0.091 


0.082/0.091 


0.081/0.091 


0.079/0.089 


0.075/0.082 




0.070/0.080 


0.079/0.089 


0.081/0.091 


0.082/0.091 


0.081/0.091 


k 4 


0.23/0.32 


0.21/0.30 


0.081/0.091 


6.36/6.93 


310/336 




310/336 


6.36/6.93 


0.081/0.091 


0.21/0.30 


0.23/0.32 


a 


0.10/0.10 


0.095/0.17 


0.081/0.091 


0.080/0.089 


0.083/0.088 


P 


0.076/0.088 


0.091/0.096 


0.081/0.091 


0.075/0.15 


0.074/0.19 


All 


0.16/0.97 


0.081/0.69 


0.081/0.091 


0.060/0.091 


0.054/0.081 


The nominal parameter values (k^ 


= 1000, /c_! = 1100, k 2 = 


1000, k- 2 = 1200, k 3 


= 1000, k. 3 = 7000, k 4 = 


1000, /c_ 4 = 1 100, a = 4.2, and p = 


1) are modified by a 



multiplicative factor, and the maximal (for any state variable) time average/infinity norm of the relative difference between the original and the reduced model is 
presented. Note that due to potential numerical inaccuracies only concentrations larger than 10" 6 are considered in the analysis. 
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state variables and not the number of parameters, which 
are reduced as a side-effect of the QSS assumptions. 

The model for glucose transport was first reduced with 
an approach in which the QSS equations and conservation 
relations were directly substituted into the remaining 
ODEs. The results of this approach are not satisfactory 
since the reduced model gives predictions that are differ- 
ent from the original model for most state variables. Our 
method was then applied to the same model both together 
with the assumption of equal concentration of the trans- 
porter in the inner and outer regions of the cell membrane 
used in [37], and without any additional assumptions. The 
application of our method together with the assumptions 
used in [37] results in a model with three state variables. 
The state dynamics is significantly better preserved than 
with the first approach, although still not satisfactory. We 
then decided to reduce the original model without the 
assumption regarding the localization of the transporter, 
with two different definitions of modified lumped state 
variables. While one of these definitions results in a 
reduced model with four state variables and gave rather 
accurate predictions, the other choice reduces the number 
of state variables to five and gives an excellent description 
of the state dynamics. It is therefore apparent that there is 
a tradeoff between accuracy and the number of state vari- 
ables in the reduction process. The glucose transport 
model corresponds to the first part of glycolysis, in which 
glucose is transported into the cell. We therefore propose 
that it might be rewarding to carefully re-investigate the 
assumptions underlying the reaction rate equations in 
complete models of glycolysis (see [43] for one example). 

We have also observed a few issues regarding the 
implementation of the method. The symbolic inversion 
of the matrix that is necessary to compute the dynamics 
of the modified lumped state variables may be expen- 
sive. However, this is typically only a practical limitation 
for large matrices, which result from large clusters of 
fast state variables. In our experience large clusters of 
fast state variables are relatively rare also in large bio- 
chemical models. Another option, if it is not practically 
feasible to invert the symbolic matrix, is to solve the sys- 
tem of linear equations in Eq. (31) numerically. We also 
observed that the symbolic right-hand side of the result- 
ing differential equations may be long. However, these 
are usually not practical limitations for the applicability 
of the method, e.g., the simulations of all examples in 
this paper are very fast on a modern computer. Avail- 
able methods to reduce the analytic reaction rate 
expressions include sensitivity analysis w.r.t. state vari- 
ables and parameters, and the method proposed in [17]. 

There is still no consensus method for automatic identi- 
fication of state variables in QSS, although criteria for the 
detection of state variables in QSS have been proposed, 
for example in [30] . A simple approach is to simulate the 



original model and investigate for which state variables the 
corresponding in- and outflow reaction rates are approxi- 
mately equal. State variables for which this condition 
holds are then considered to be in QSS. Note that for the 
models in this paper it was already clear from the bio- 
chemical understanding of the corresponding systems 
which of the state variables that could be considered fast 
(see [37] for the glucose transport example). However, an 
appropriate general criterion for automatic identification 
of state variables in QSS is still lacking. 

Although the theory presented in this paper constitutes 
a great leap forward for construction of zoomable mod- 
els, more research is required to make the method fully 
automatic. An important challenge is to define a mean- 
ingful measure for the similarity between the hierarchical 
model layers (degrees of zooming). Another interesting, 
although trivial, observation that deserves further atten- 
tion is that QSS assumptions typically do not hold in the 
whole parameter space. Although the reduced models in 
this paper appear to be robust to varying parameter 
values it may not be the case in general. It may therefore 
be revealing to compare the original model and the 
reduced model to characterize the parameter space 
regions in which the QSS-assumptions are valid. 

Conclusions 

We have presented a novel method for reduction of bio- 
chemical models that is compatible with the concept of 
zooming. Zooming allows the modeler to operate on 
different levels of model granularity, and enables a direct 
interpretation of how modifications to the model on one 
level affect the same model on other levels in the hierar- 
chy. The proposed method is based on the application 
of proper lumping in combination with the identifica- 
tion of linear relations in nonlinear equations. 

The method was applied to two example models. The 
first model is small and commonly occurring as a part 
of larger biochemical models. The second example is a 
model for glucose transport in baker's yeast, which con- 
stitutes the starting point for glycolysis. Both models 
could be significantly reduced with the proposed 
method, and the resulting state variables could be back- 
translated to the original state variables. The method 
that is presented in this paper constitutes an extension 
of the method that was previously developed for linear 
biochemical models to its nonlinear counterpart. Since 
most models in the systems biology community are in 
fact nonlinear, our method constitutes an important 
step towards zoomable biochemical models. 

A Appendix 

A.I Appendix 1 

The ordinary differential equations for the enzyme 
kinetics model take the form 
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where r x = k^SE - k^C s , r 2 = k 2 C s , and r 3 = k 3 C P - k. 
3PE. 

The parameters are set to values that satisfy the 
assumptions of dominating and insignificant reaction 
terms, with k x = 1000, k_ x = 2000, k 2 = 1, k 3 = 3000, and 
/c 3 = 1000, together with the initial conditions 5(0) = E 
(0) = 1 and P(0) = C s (0) = C P (0) = 0. This gives the 
parameter values Mi = 0.5, M 3 = 3 and L E = 1 in the 
reduced model. Initial conditions can in general be 
obtained from a short simulation of the original model, 
until the fast state variables reach QSS, but in this case 
Eq. (19) gives an analytic expression of 5(0) (P(0) = 0) 



1 



5(0) = -[L s -L E - M- 1 + ^(L s + L E + M~ l Y - 4L S L E ) 
= V3 - 1, 

where M^ 1 = K\ , S = A, L s = L lf and L E = L 2 . 



A.2 Appendix 2 

The ordinary differential equations for the glucose 
transport model take the form 

dx e Gl 

= ~kiX e E x e Glc + k-ix? E _ Glc , 



dx\ 



dt 



dt 

dx' 



-jj^ - —k 2 X E X l Glc + k- 2 X E -Glc' 



- k4X l E X l G6 p k-4%E-G6P' 

= k3X l E _ Glc X l G6 p — k-3X E _ Glc _ G6p , 



dt 



dt 



dt 



dt 
dt 



- k3 x E-Glc X G6P + k-3 X E-Glc-G6P 
k4X E X G6p + k-4X E _ G6p , 

- = oi(x l E _ Glc — *£_ G j c ) + kiX e E x e Glc — k-\X E _ Glc , 

- = Oi(x E _ Gk — X E _ Gk ) — k3X l E _ Glc X l G6 p + 
k-3 X E-Glc-G6P + ^ X E X Glc ~ ^-2 X E-Glc' 

= fi{x E — x e E ) — kiX e E x e Glc + k-iX e E _ Glc , 

= fi{x E — X l E ) — k4X l E X l G6 p + k-4X l E _ G6p — 



k 2 x l E x Glc + k- 2 x x E _ Glc , 



which was introduced in [27]. 
A.3 Appendix 3 

In this section we investigate an alternative (naive) 
approach to reduce the glucose transport model. The 
first step is to identify state variables for which the QSS 
assumption holds, and the mass conservation relations 
in the model. In the second step of this approach we 
then substitute the corresponding system of equations 
into the ODEs corresponding to slow state variables. 

Now consider the model for glucose transport in 
yeast. We assume that the state variables x? Glc , x l Glc , 

x e-G6P an< ^ x E-Glc-G6P are m Q^S> which gives Eqs. 
(43)-(46). Note that Eqs. (45)-(46) indirectly imply that 
x l G6P is in steady state. The substitution of Eqs. (43) - 
(46) into the ODEs of the original model gives 



Glc 



dt 



a i X E-Glc X E- 



Glc 



dx l 



E-Glc 



dt 



a{xf 



E-Glc 



x E-GlcJ 



dt 



4)- 



(61) 



(62) 



(63) 



(64) 



Note that the state variables x \-q\ c and x l E _ Glc are 
decoupled from the state variables x e E and x l E in Eqs. 
(61)-(64). 

There are three molecules (moieties) whose mass is con- 
served in the model as a whole, i.e., Glc, G6P, and E. How- 
ever, we can not substitute any of the conservation 
relations into the remaining ODEs without re-introducing 
state variables that were already eliminated. So the final 
reduced model takes the form of Eqs. (61) - (64). However, 

the sum of the state variables x e-ci c an d x e-gIc> an< ^ x e 
and x l E is conserved in the reduced model, which makes it 
possible to reduce the model to two state variables. 

Unfortunately, due to the form of the ODEs and the 
initial conditions of the state variables in the reduced 

model, the state variables X E -Glc an d x e-gIc remam equal 
to zero at all times, and only the state variables x e E and 
x l E take non-zero values. We therefore decided to simu- 
late the original model for a short time until the fast state 
variables reach QSS, and to use the final state variable 
values as initial conditions in the reduced model. 
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The solution to the equation system defined by Eqs. 
(43)-(46) and Eq. (49) is 



X Glc ~ K l ■ 



X 



E-Glc 



^Glc 



X 



,i 

E-Glc 



K3K4 



X G6P 



X E-G6P 



X E-Glc-G6P 



K4X l g_Qi c + K^x 1 ^ + K3X4 
K^x l E 

^4^7_ GZc ^3 X E ^3^4 



^G6P/ 
^G6P/ 
^G6P/ 



which can be used for back-translation of the state 
variables of the reduced model to those of the original 
model 

The predictions of the state variables of the original 
model, resulting from simulations of the original model 
and the reduced model with the parameter values set as 
in [37], is not satisfactory for any other state variable 
than G6P, which remains approximately constant over 
time. Implementations of the original model and the 
reduced model in SBtoolbox2 for MATLAB [36] are 
included in Additional files. 

A.4 Appendix 4 

In this section we apply our method to the glucose 
transport model, and following [37] we will assume that 
the concentrations of the transporter are constant in the 
inner and outer regions of the cellular membrane. With 
this assumption the distribution among the transporter 
state variables of the original model, which constitute 
the lumped state variables L E , is uniquely defined. 

The first step of the method is to identify the apparent 
conservation relations in the model. We note that G6P 
and the transporter E are conserved, and apparent con- 
served glucose (see Definition 1) in the inner and outer 
regions of the membrane, respectively. The four appar- 
ent conservation relations take the form 



: MX ■■ 



Lac 2 

LG6P 

\Le / 

(\ 0 0 0 0 1 0 0 0\ 
0 10 10 0 10 0 
0 0 1110 0 0 0 

\0 01101111/ 



x, 



where 



X £ X e) T 



In the second step of the method we define the modi- 
fied lumped state variables. We decide to keep L E in the 
reduced model since it corresponds to an exact conser- 
vation, and therefore results in the largest reduction 
possible (note that the exact conservation relations, L E 
and L G6P , can not simultaneously be used since the 
lumping would then not be proper). The modified 
lumped state variables take the form 



V L E J 



-- M m x -- 



/100000000\ 
010000000 
000010000 

Voo 1 1 0 1 1 1 1/ 

We note that Eqs. (43)-(47), and (50) are all linear w.r.t. 
the state variables that constitute state L E , so the require- 
ment for the existence of at least 6 (n m ) linear relations is 
satisfied, which enables back-translation in step three of 
the method. 

In the third step of the method we derive the fraction 
parameters for the lumped state variable L E . Eqs. (43)- 
(47), and (50) form an equation system, corresponding 
to Eq. (28) 



(0 0 1 
000 
1 0 0 
0 10- 
00a 

M 1 1 



1 



0 \ 



" A G6P 
0 



1 / 



*E-G6P 

i 

E-Glc-G6P 
X E-Glc 
X E-Glc 



4, 



/ 



0 
0 
0 
0 



(65) 



where of Glc 



fe_i X Glc > 



x Glc ~ k- 2 Glc> ^G6P 

X G6P = fe^4 x G6P* ^ ne s °l ut i° n to Eq. (65) is given by Eq. 
(29) 



/ 4. 

X E-Glc-G6P 
X E-Glc 
X E-Glc 
X% 



\ ( fi-G6P \ 
^E-Glc-G6P 
V E-Glc 
lE-Glc 

\ Ve 
( X G6P(P + aX ilJ 



) 



X G6P X Glc(P + aXe GlJ 
S 

x e Glc {/3 + ax l Glc ){/3 + ax e Glc ) 



(P + ax Glc )(/3 + ax Glc ) 



(66) 
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where 

£ = (P X Glc + 2^fc a 4;fc + 2 ^ + a 4fc + aX Glc 



X Glc X G6P ' 



+ 0(X l Glc X l G6p xf Glc + f3x l Glc + fix l G6p + G^gZc X G6p) 



(67) 




and where the fraction parameters were calculated 
with Eq. (30). We note that the fraction parameters are 
functions of xf Gkf x l Glc , and x l G6P , which are state vari- 
ables both in the original- and in the reduced model. In 
the fourth step of the method we derive differential 
equations for the modified lumped state variables. The 
ODE for the fourth state is l m = L E = 0, which is 
replace by a constant. The ODEs for the other states are 



J 12 7l3 

1 + J22 J23 

J32 1 + 733 , 



where l mv3 denotes the first three states variables in 
l m . Note that there are three state variables in the 
reduced model, which is the same number as for the 
reduced model in [37]. 

In the fifth step of our method we compare predic- 
tions of the original state variables between the original 
and reduced models, where L E is back-translated with 
the fraction parameters defined in Eqs. (66) -(67). 

The simulation results are clearly more accurate than 
with the approach in Appendix A.3, although still not 
satifying. We refer to Additional files for implementa- 
tions of the original and reduced models in SBtoolbox2 
for MATLAB [36]. 

A.5 Appendix 5 

In this section we apply our method to the glucose 
transport model, but with an alternative definition of 
the modified lumped state variables. We do not use the 
assumption of constant regional concentrations of the 
transporter (Eq. (50)). 

In the first step of the method we note that the appar- 
ent conservations are given by Eq. (51). 

In the second step of our method we decide to keep 
state variable Le 2 , instead of L G6P , in the reduced model. 
This leads to the following definition of the modified 
lumped state variables 



= M m x = 



/ X Glc \ 



A Glc 
X G6P 

\Le 2 J 

/100000000\ 
00000 10 10 
010000000 
000010000 

\00 1 1 00 1 0 1/ 



(68) 



Where 

X = ( X G/c X Glc X E-G6P X E-Glc-G6P X G6P ^E-Glc X E-Glc X £ X eY > 

and we note that x e Gk , x l Gk , and x l G6P are state variables 
both in the original and reduced models. Also note that 
the requirement of at least n m equations, that are linear 
w.r.t. the original state variables and linearly indepen- 
dent, is satified for each of the modified lumped state 
variables by Eqs. (43)-(49). 

In the third step of the method we calculate fraction 
parameters for the modified lumped state variables L El 
and L El , which correspond to more than one of the ori- 
ginal state variables. The fraction parameters for state 
variable are given by Eq. (53). We can now use Eqs. 
(44) -(46) and Eq. (68) to form an equation system corre- 
sponding to Eq. (28) 



A(l m ,p)x 



K 2 
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0 -X* 



0 0 \ 

K 4 0 
, 0 K 3 
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/ 0 \ 
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VW 
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= b k (lm,p), 



\ 



A E-Glc 



X E-G6P 
\ X E-Glc-G6P , 



with the solution given by Eq. (29) 
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(69) 



Le 2 > 



/ 



where 



X G6P 



(K 4 x[ 



Qi c + + K^K^x l Glc + K2K3K4 



and the fraction parameters were calculated with Eq. 
(30). 

The fourth step of the method is to derive ODEs for 
the modified lumped state variables. Since the apparent 
conservations are separated into two disjoint clusters of 
fast state variables, we can treat the model for the inner 
and outer regions of the membrane separately. The rate 
equations for the outer region are given by Eqs. (56)- 
(57). Eq. (31) gives us the ODEs of the modified lumped 
state variables in the inner region 
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where / 3:5 and \ my5 are the last three state variables of 
/ and / m , respectively, and 



'3:5 



LG6P J 

v Le 2 / 



-ot{x\ 



E-Glc ~ 
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x E-Glc 



(71) 



a ( X E-Glc X E-GlJ P( X E X e) , 



In the fifth step of the method we simulate the 
reduced model with Eqs. (56)-(57) and (70)-(71) and we 
then use Eqs. (53) and (69) for back-translation of the 
state variables. A comparison between the original 
model and the reduced model, w.r.t. the state variables 
of the original model, is presented in Figure 5 and Fig- 
ure 6. The agreement between the models is very good. 
We refer to Additional files for implementations of the 
original and reduced models in SBtoolbox2 for 
MATLAB [36]. 

Additional files 

The original and reduced versions of the models pre- 
sented in this paper, and scripts for simulation and 
comparison between the original and reduced versions 
of the models. Note that the systems biology toolbox for 
MATLAB [36] and the symbolic math toolbox for 
MATLAB must be installed on the system for simula- 
tion of the attached models. 

Additional material 



Additional file 1: Model 1. The original enzyme kinetics model. 

Additional file 2: Model 2. The reduced enzyme kinetics model. 

Additional file 3: Script 1. Script for comparison between the original 
enzyme kinetics model and the reduced model. 

Additional file 4: Model 3. The original glucose transport model. 

Additional file 5: Model 4. The reduced glucose transport model with 
the alternative (naive) approach. 

Additional file 6: Script 2. Script for comparison between the original 
glucose transport model and the reduced model with the alternative 
(naive) approach. 

Additional file 7: Model 5. The reduced glucose transport model with 
our method and the assumption of constant concentrations of the 
transporter in the inner and outer regions of the cellular membrane. 

Additional file 8: Script 3. Script for comparison between the original 
glucose transport model and the reduced model with our method and 



the assumption of constant concentrations of the transporter in the 
inner and outer regions of the cellular membrane. 

Additional file 9: Model 6. The reduced glucose transport model with 
four state variables with our method. 

Additional file 10: Script 4. Script for comparison between the original 
glucose transport model and the reduced model with four state 
variables with our method. 

Additional file 11: Model 7. The reduced glucose transport model with 
five state variables with our method. 

Additional file 12: Script 5. Script for comparison between the original 
glucose transport model and the reduced model with five state variables 
with our method. 
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